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ABSTRACT 

Orbiting disks may exhibit bends due to a misalignment between the angular 
momentum of the inner and outer regions of the disk. We begin a systematic 
simulational inquiry into the physics of warped disks with the simplest case: the 
relaxation of an unforced warp under pure fluid dynamics, i.e. with no inter- 
nal stresses other than Reynolds stress. We focus on the nonlinear regime in 
which the bend rate is large compared to the disk aspect ratio. When warps are 
nonlinear, strong radial pressure gradients drive transonic radial motions along 
the disk's top and bottom surfaces that efficiently mix angular momentum. The 
resulting nonlinear decay rate of the warp increases with the warp rate and the 
warp width, but, at least in the parameter regime studied here, is independent of 
the sound speed. The characteristic magnitude of the associated angular momen- 
tum fluxes likewise increases with both the local warp rate and the radial range 
over which the warp extends; it also increases with increasing sound speed, but 
more slowly than linearly. The angular momentum fluxes respond to the warp 
rate after a delay that scales with the square-root of the time for sound waves 
to cross the radial extent of the warp. These behaviors are at variance with a 
number of the assumptions commonly used in analytic models to describe linear 
warp dynamics. 

Subject headings: accretion, accretion disks - hydrodynamics 



1. Introduction 



Astrophysical manifestations of warped accretion disks can occur at a myriad of scales 
in the Universe, requiring no more than a misalignment of angular momentum between the 
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inner and outer regions of the disk. The exact cause of this misalignment may vary. One of 
the first analyses of misaligned fluid disk dynamics was in the context of a Kerr black hole 
whose spin is oblique to the rotation o f the disk; the Lense-Thirring effect then forces disk 
precession faardeen fc Pettersonlfl975T ). Alternatively, a disk orbiting a proto-stellar binary 
at a tilt to the binary o rbital plane may develop a warp due to the quadrupol ar term in the 
gravitational potential ( jPapaloizou fc Terquemlll995t iLubow fc Ogilvid |2000| ) . Even in the 
absence of an external torque, the orbital plane of the mass supply may change over time, 
e.g. in the case of proto-stars, AGN, and the specific cases of SS433 and Her X-l. Absent 
external effects entire ly, disks may undergo self-induced warping due to their own radiation 
output (jPringlelll996[ ). The presence of these warps is likely to have important consequences 
for the internal dynamics of the disk as well as observational significance. 

By definition, if a disk is warped, the orientation of orbital angular momentum in the disk 
changes with radius. Thus, in the end, whether warps grow, decay, or mutate, their evolution 
is very largely dependent on how angular momentum in different directions is delivered to 
the disk (e.g., via the Lense-Thirring mechanism) and then moves through it (i.e., by the 
action of internal stresses). Torques due to external causes can often be evaluated directly, 
but internal stresses are a much harder problem. Although th ey are now known to be due to 
MHD turbulence stirred by the magneto-rotational instability (IBalbus & Hawleylll998l ). even 
in the case of flat disks it is difficult to quantify these stresses in detail without large-scale 
numerical simulations. Little is known about the character of this sort of MHD turbulence 
in warped disks; it is quite possible that its nature changes in significant ways. At the 
very least, in the fully three-dimensional context of warped disks, the number of interesting 
components in the stress tensor grows from the single one relevant to flat disks (the r-0 
component) to all six independent components, particularly the other two off-diagonal ones. 
The problem of describing the internal stresses is further complicated by the fact that disk 
bending induces ra dial motions because the be nd misaligns the vertical pressure profiles of 
neighboring rings ( jPapaloizou fc Pringldll983l ). These contribute to the r-z component of 
the Reynolds stress precisely because bending means that the disk midplane changes as a 
function of radius. However, it is not easy to estimate the magnitude of this contribution to 
the stress because these radial motions can be limited by a variety of mechanisms, such as 
the work associated with fluid compression, changes in the pressure profiles over time, and 
the MHD stresses. 



Because of these difficulties, the overwhelming m ajority of work in this field h as instead 
described internal stresses in terms of the ansatz of IShakura fc Sunyaevl (119731 ). Analyz- 
ing the radial structure of time-steady, unwarped, thin disks, they pointed out that the 
internal stress tensor has only one interesting element, the one responsible for conveying 
the component of angular momentum parallel to the disk's mean angular momentum in 
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the radial direction (i.e., the r-cf) component). Arguing on the basis of dimensional analy- 
sis, they suggested that its magnitude should be ap, where p is the local pressure in the 
disk and a is a parameter thought to be of order or somewhat less than unity. In subse- 
quent treatments, this stress has often been modeled as being due to a phenomenological 
isotropic anomalous viscosity, s o tha t it is proportional to the local fluid shear rate. For 
example, iPapaloizou fc Pringld ( 119831 ) developed a one-dimensional formalism for describing 
the evolution of linear disk warps in which they supposed that the internal stress for all three 
off- diagonal components was due to a phenomenological viscosity proportional to the local 
pressure and the appropriate fluid shear. 



It is convenient 



;o di vide this earlier work according to the severity of the warp. 



Nelson fc Papaloizoul (119991 ) pointed out that the warp becomes nonlinear when the bending 
rate dO/dhir > H/r, where 9 denotes the angle of the local orbital normal from a reference 
direction and H is the disk scale height. Bends sharp enough to make d8/d\nr > H/r (as 
we will discuss later) induce transonic radial motions, so that the nonlinearity is in the sense 
that the perturbed velocities in the disk are comparable to or greater than the sound speed. 
Because most accretion disks are expected to be quite thin, any substantial tilt between the 
disk's orientation at large radius and its orientation at smaller radius would be stretched 
over many decades of radius unless the warp is nonlinear. For this reason, nonlinear warps 
are an important topic. 

Nonetheless, equations are always much more tractable in the linear regime, so a large 
part of the analytic effort devoted to warps in disks has been restricted to the limit of 
gentler warps. This linear regime is generally further subdivided into two limiting cases: the 
diffusive limit, in which a > H/r; and its opposite, the bending wave limit. In both cases, 
the warp still drives radial motions, but they are subsonic. They are, however, generally 
treated differently. In the diffusive case, the associated Reynolds stress is seen as driving 
warp relaxation, but the isotropic anomalous viscosity limits its amplitude. The end-result 
is th at warp relaxa tion is modeled by a diffusion equation whose diffusion coefficient cti oc 
oT x (IPringld Il992l ) . In the long- wavelength bending wave limit, rather than causing warp 
relaxation, the radial motions are thought to create a circulatory motion in the poloidal 
plane which drives a propagating wave with speed ~ c s q/2, where c s q is the isothermal sound 
speed. The rol e of isotropic anomalous viscosity is then restricted to a slow damping of the 



bending wave (IPapaloizou fc Lin 



1995 



Demianski fc Ivanov 1997; Ivanov fc Illarionov 1997 



Lubow fc Qgilvid l2000l ; iLubow et al.ll2002l ). In fact, it is this estimate of the bending wave 
damping rate that underlies the d istinction between these two regimes according to the ratio 
a/(H/r) JPapaloizou & Linlll995h . 



By separating the dynamics in terms of velocity scale (orbital vs. transonic radial 
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motions vs. net inflo w), performing an asymptotic expansion in H/r, and then integrating 
over spherical shells, lOgilvid (119991 ) was able to transform a nonlinear description of the 
fluid motions into a 1-d equation for the radial evolution of the disk's local orientation. In 
the course of this transformation, he derived expressions for the «2 parameter. However, 
extending this formalism into the bending wave regime has proved difficult. One expression 
of this difficulty is the fact that no solution for a 2 can be found when a = 0, the e picyclic 
frequ ency k is less than the orbital frequency Q, and the bending rate d6 / lnr > 0.2 (jOgilvie 
19991 ). Unfortunately, this range of parameters is, in fact, generic for nonlinear bending waves 
in the absence of internal stresses because hydrostatic equilibrium with non-zero pressure 
always makes k < tt. 

Given the obstacles to analytic work in nonlinear fluid dynamics, many have also 
turned to numerical simulation. The overwhelming majority of simulations performed thus 
far has used the smoothed particle hydrodynamics (SPH) method including an isotropic 



anomalous viscosity. These simulations 
warps produced by a binary potential ( 



rave explored the dy namics of warp evolution for 



Larwood et al.l Il996l ). warps produced through a 



Lense-Thirring torque dNelson fc Papaloizou 2000 ), and the relaxation of a n unforced warp 
jNelson fc Papaloizoul Il999l : lLodato fc Prinrie! 120071 : lLodato fc Pricel boioh . Their results 
are generally consistent with the linear theory (likewise including an isotropic anomalous 
viscosity) in both the diffusive and bending wave regimes. Some nonlinear effects have also 
been explored. 

However, there are still a number of ways in which our understanding of this phe- 
nomenon remains unsatisfactory. In particular, nearly all of the extant results depend on a 
model for the disk's internal stresses that does not correspond directly to the actual physical 
mechanism. Unlike the a-model of stress, MHD stresses are not necessarily linked to the 
local pressure. Unlike the stress due to an isotropic shear viscosity, it is the time-derivative 
of the Maxwell stress that is related to the shear, not the stress itself. Moreover, the na- 
ture of the relationship of stress to shear is different from that of shear viscosity, is highly 
anisotropic, incorporates proportionality factors quadratic in the magnetic field itself, and 
depends on gradients of the magnetic field in addition to gradients of the fluid velocity. Thus, 
it is unclear to what degree either an isotropic phenomenological viscosity or a diffusivity of 
the a 2 variety mimics the effects of these turbulent Maxwell stresses. 

In order to answer these questions, we are embarking on a program in which we will 
approach the problem of internal stresses in warped disks from a different point of view. We 
will consider only effects derivable directly from hydrodynamics and magnetohydrodynamics. 
To be able to identify which mechanisms are responsible for which effects, we will introduce 
them one by one. In this paper, we wish to study the effects of the fluid motions induced 
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by the warp in isolation — without MHD and also without any phenomenological viscosity 
to artificially limit them. Thus, this work might be thought of as studying the nonlinear 
bending wave regime because a = 0; put another way, numerical techniques permit study of 
exactly the interesting parameter regime in which the extant nonlinear analytic theory fails. 
With the data from numerical hydrodynamics simulations, we can measure the Reynolds 
stresses driven by warps and explore their dependence on initial warp strength and physical 
parameters of the disk such as sound speed. We will introduce MHD effects in a later paper 
now in preparation. 

The plan of the paper is as follows: Section 2 describes our numerical method and the 
properties of the warped disks we simulate; Section 3 presents an overview of how disks with 
nonlinear warps relax when only hydrodynamic mechanisms are active; Section 4 analyzes 
the internal stresses in greater detail; and Section 5 summarizes our results and places them 
in the context of previous work. In Appendix A we describe how we tested our methods for 
numerical dissipation. 



2. Warped Disk Model 
2.1. Equations, initial and boundary conditions, and grid definition 

We study the hydrodynamical mechani sms involved in warp re laxation through simu- 



lations conducted using the Zeus algorithm ( jStone fc Norman! Il992l ). We use an implemen- 



tation of Zeus that solves the three-dimensional equations of hydrodynamics in coordinates 
that have a diagonal three-metric, gu- In this formulation, the equation of mass conservation 
is 

! + -La i( VW) = o ID 

where p is the density, 7 is the determinant of the three-metric, and u l is the contravariant 
velocity component. The j component of momentum density is defined in terms of the 
covariant momentum density, pwj, so that the evolution equation for the momentum density 
is written as 

^ + -±= d t {Vipwju*) + d j (P + Qjj) - ^-d j9kk - pd s * = (2) 

where P is the pressure, and <3> is the gravitational potential. The derivative of the metric 
components accounts for the inertial forces associated with a choice of coordinate system 
(e.g., centrifugal and Coriolis forces). The symbol Q denotes the stress tensor associated with 
the artificial bulk viscosity required for proper treatment of shocks; it is diagonal, and Qjj is 
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the element associated with direction j. We use an ideal gas equation of state, P — (Y — 1) e, 
where e = pe is the internal energy and e is the specific internal energy. The internal energy 
equation is 



The magnitude of the artifici al viscous stress Q a is ba sed on the convergence of the physical 



velocity components, v 1 as in IStone &; Norman! ( 119921 ) 



Q t = {n s Axf p {d t vf , (4) 

when diV % < (no sum over i); Q = elsewhere. The constant n s determines the number 
of zones over which the artificial viscosity spreads out a shock. In Zeus applications this 
is typically n s = 2, which we use here. This numerical viscosity thermalizes kinetic energy 
lost in the shocks by increasing the entropy of the gas; by following entropy generation 
we can thus determine when shocks are an important dynamical element of an evolution. 
The artificial viscosity also reduces zone-to-zone oscillations that would be excited by such 
shocks and improves code stability. By design, the artificial viscosity has negligible (or zero) 
magnitude except in strongly compressive regions, nor does it create any shear viscosity. 

Because we are studying warped disks, there is no way to construct a single system 
of coordinates in which the orbital velocity consistently points parallel to a grid coordi- 
nate. Consequently, greater numerical diffusion can be expected than in simulations of flat 
disks. On the basis of the tests we discuss at greater length in Appendix |Aj we decided 
that spherical coordinates were the best choice for conducting our simulations. In spher- 
ical geometry, even when disks are tilted by as much as 30° relative to the coordinates, 
numerical artifacts are held to a very low level when the gridscale gives a resolution of at 
least ~ 8 Zones Per vertical scale Height (ZPH). For all our simulations, the computational 
domain is (r,9,4>) e [1, 19] x [0.05, 0.95]7r x [0, 2tt]. The unit of length is arbitrary because 
Newtonian gravity has no characteristic scale. The size of the radial cells is logarithmically 
graded, with Ar increasing outwards, and the azimuthal variable, 0, is uniformly spaced. 
Both Ar/r and A<p are set equal to AO so that the spatial resolution is isotropic. For 
the HW series of simulations (defined below), the polar angle, 0, is also uniformly spaced; 
however the increased computational cost would make this uniform spacing infeasible for 
the thinner T W simulation . Inste ad, the polar angle is spaced using a polynomial spacing 



(Equation 6 of iNoble et al.l ( 120101 ) . with £ = 0.7 and n = 11) to ensure adequate resolution 
near the midplane while accepting lower resolutions at higher altitudes. For this simulation, 
A(j) = 2A6(6 = 7r/2) = 2 Ar/r. The resolutions used for the simulations presented here have 
(N r , Nq, JV ) = (128, 128, 192) and (N r , N e , N+) = (256, 128, 240) grid cells for the HW and 
TW runs, respectively. These simulations are run using outflow boundary conditions in the 
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radial and polar coordinates and a periodic boundary condition in the azimuthal direction. 
The lack of a magnetic field in these simulations implies that there will be no consistent 
stress to drive inward flow, and thus there is little interaction between the disk and the 
boundary. 

To avoid both dynamical transients associated with the disk boundaries and numerical 
effects associated with the boundary conditions of the simul ation domain, we use as our initial 



condition a state derived from an exactly hydrostatic torus (lHawleyll2000l ). Warps, of course, 
introduce non-hydrostatic features, so our simulations do not begin from an equilibrium. Our 
initial conditions are constructed in two steps: first we lay down a density, pressure, and 
velocity distribution consistent with a flat hydrostatic torus; then we apply a systematic 
warp. The result is that the angle Bt between the local orbital normal (the direction of 
the total angular momentum at a particular spherical radius r) and the z-axis is a specified 
function of r. We define the x-axis so that the tilt is in the x-z plane. 

A hydrostatic torus is fully determined by five parameters: the shear parameter q de- 
termined by the non-Keplerian orbital velocity profile (Q oc R~ q ); the adiabatic index T; 
the radius of the inner edge, i?j n ; the radius of maximum pressure, Rm] and the maximum 
density, p(R = Rm, z = 0) = pu- The units of density are arbitrary because the gas is not 
self-gravitating. Because of its axisymmetry, the properties of such a torus are functions of 
only R and z (cylindrical radius and height from the midplane), i.e. p(R,z), e(R,z), and 
v^(R) = RQ(p. These dependences can easily be translated to the spherical coordinates of 
the simulation by the usual relations. Its midplane is, of course, located at z = 0. 

In the work presented here we use two different initial tori. For both, T = 5/3, R{ n = 2, 
Rm = 4, and pm = 100. However, they differ in their values of q. For our fiducial model, 
designated (H)ydrodynamic (W)arped disk (HW), q = 1.6; in our alternate model, called 
(T)hin (W)arped disk (TW), q = 1.52. As a result, they differ in their typical sound speeds 
and aspect ratios. In the HW torus, H/r w 1/4 throughout the body of the disk, while 
the TW torus is more nearly Keplerian and therefore also has less pressure support, so that 
H/r « 1/10. 

To warp the torus, we first choose a function 9t{t). In order to isolate the warp in the 
middle of the disk and ensure continuous variation, we use functions of the form 




r <r c - L w , 

\r-r c \<2L w (5) 
r > r c + L w 



where 9p is the angle defining the orbital inclination of the outer disk, and the coefficients 
A and B are chosen so that ^(r) is continuous. The midpoint of the bend r c is set to 
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coincide with the radius of the spherical shell bounding half the total angular momentum of 
the system. The central radius r c = 9 for the HW series and r c = 8.3 for the TW simulation. 

With this function specified, we define a warped cylindrical coordinate system (R', </>', z') 
related to the original cylindrical coordinate system through rotating a Cartesian coordinate 
scheme by 9t{t) about the ?/-axis and then using the usual definition (i.e., R' = (x' 2 + y' 2 ) 1 ^ 2 , 
(f)' = taxi^ 1 (y' / x')) . The data defining the original hydrostatic torus are then remapped onto 
this new coordinate system by identifying points: p(R', z') = p(R, z), etc. The magnitude of 
the orbital velocity is preserved in the midplane, but its direction is, of course, made parallel 
to <p' . Off the midplane, the warped disk is constrained to orbit on cylinders. This warped 
coordinate system is used only for construction of the initial condition; the simulation itself 
is conducted in spherical coordinates. 



2.2. Quantifying the warp: t/j 

To complete the definition of a warped disk, we must choose the parameters A, B, and 
dp. We do so by characterizing the warp in terms of the dimensionless quantity 

dO T {r) 

m = -^r. (6) 

The functional form we have chosen results in i/j(r) being piecewise constant. Note, however, 
that our use oiip in this context differs slightly from the standard definition for if) (= r\d£/dr\, 
with i a unit vector pointing in the direction of the local angular momentum), in that 
Equation [6] is a signed quantity and measures only rotation about the y— axis. In the case 
of negligible £ y , our quantity agrees in magnitude with the standard definition. 



As first remarked by lNelson fc Papaloizoul (119991 ) , the linear regime in a disk with H/ r <C 



1 is defined by ip <C H/r, suggesting that ip = ip/ (H/r) is also an important measure of disk 
warp. To see why this is so, we will use analytic estimates here and reinforce these points 
quantitatively when we describe our numerical results. 

The fundamental reason why ip > 1 marks the onset of nonlinearity is that this is the 
criterion for the vertical displacement across a radial separation ~ r to be larger than a scale 
height H . When this is so, the warp creates a radial pressure contrast at fixed height from 
the midplane that is order unity. As a result, fluid forces in the radial direction are no longer 
small perturbations to the flow. 

In addition, when ip > 1 across a radial extent ~ r or more, the elongation of the 
pressure contours in the warped region becomes significantly offset from the radial direction 
(see Fig. [2]). That is, the plane of the pressure distribution does not coincide with the local 
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orbital plane. In a thin disk, the local pressure gradient lies very close to the normal to 
the equatorial plane, but the pressure gradient due to a warp tilts the pressure gradient 
oblique to the orbital normal. When ip > 1, this obliquity becomes sizable. Consequently, 
the matter cannot remain hydrostatic even in the rotating frame. Moreover, because the 
pressure gradient (whose natural length scale is H) is not balanced by gravity, the fluid 
velocity it induces can be large. Integrating the acceleration due to the pressure gradient 
over an orbital period leads to transonic radial velocities: 

5v r ~ P/(pHQ) ~ c„ (7) 

where c s is the fluid sound speed. 

For these reasons, in this paper we will be primarily interested in initial warps for which 
r/S > 1 where the disk bends. Although H/r is not exactly constant throughout our disk, it 
varies slowly enough that ip is nearly proportional to ip, and the warp can be characterized 
in terms of only two parameters, ip and the radial width of the bend region, 2,Lw- 

The parameters for the suite of simulations presented here are given in Table HJ where 
the simulation ID encodes the family of simulation, HW and TW for H/r m 1/4 and 1/10 
respectively, and also gives the approximate value of ip in the intermediate radial range as 
well as the radial domain over which the disk is warped. We use the notation (V)ery (W)ide, 
(W)ide, and (N)arrow to refer to L w = 6,4,2 respectively. Due to the reduced computational 
cost of thicker disks, we assembled a suite of HW runs of varying warp amplitudes and radial 
spans, whereas we limited ourselves to only one TW simulation. TW is a thinner analogue 
of HW-1 . 5W: its warp amplitude and radial span are identical, but the thinner disk results in 
a larger value of ip. 

The initial radial profiles of ip for all the simulations are shown in Figure [IJ where it 
can be seen that there is little radial variation in if> in the warping region of the disk. An 
example of the warped torus model is presented in Figure [2j in which the density contours 
of the disk at = are shown superimposed on a map of the warped coordinate system 
(R ; , z') for simulation HW-1 . 5W, which we will use as a fiducial model. 

We show an example of the end-result of this procedure in Figure |2l where the initial 
structure for simulation HW1.5 is shown. 

All the simulations were run for 10 orbital periods at the midpoint of the warp, r c , and 
we use this orbital period as the unit of time. Thus, in these units the sound speed in the HW 
simulations is ~ 14 and in the TW simulation 5.6. Full 3-D datasets were recorded every 
0.1 period. When we quote averaged values below, we denote a volume- weighted average by 
(...); a density-weighted average is indicated by the notation (. . .) p . 
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Run ID 


Of 


L w 


Approximate ip 


HW-0 . 5VW 


12° 


6 


0.5 


HW-1.5W 


20° 


4 


1.5 


HW-2 . 5W 


35° 


4 


2.5 


HW-2 . 5N 


15° 


2 


2.5 


HW-5N 


35° 


2 


5 


TW 


20° 


4 


3.75 



Table 1: Simulation parameters. 




Fig. 1. — Initial radial profile of ip(r) for all simulations. 
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X-Axis 

Fig. 2. — Density contours of a slice in the x — z plane through a warped disk (color 
curves) shown relative to two coordinate systems. The axes are marked in terms of Cartesian 
coordinates (x,z), while the black lines show the grid of the warped cylindrical coordinate 
system (R', z') in this plane. Note that the (R', z') coordinates are not orthogonal in the 
warped region of the disk due to the non-zero derivative of 8t, however a local orthonormal 
basis can be defined. The lack of orthogonality does not concern us here as these coordinates 
are merely used to define our initial conditions. 
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3. Results 

3.1. Radial Pressure Gradients and Radial Streaming 

In any hydrostatic flat disk configuration, there are finite radial pressure gradients that 
balance the gradient in the effective potential due to the slightly non-Keplerian angular mo- 
mentum distribution. As expected, in our warped disk there are unbalanced radial pressure 
gradients right from the start. Figure |3] shows the residual pressure gradient, dlnP'/dr at 
radius r = 9 (the middle of the warp region) at t = in both the always-linear simulation 
HW-0.5VW and the modestly nonlinear simulation HW-1 .5W. The residual pressure gradient is 
defined as 

dlnP' _ dlnP d\nP 

dr dr dr 

where d\nP /dr is the midplane value of d In P/dr in an unwarped torus. Within 1-2 scale 
heights of the midplane, dlnP^/dr changes minimally, so it is also a very good approximation 
to the radial pressure gradient even well off the midplane. At r = 9, we find dlnP^/dr m 
—0.4. We make this adjustment to the pressure gradient in order to emphasize the portion 
of it due to the warp. 

In the linear simulation, we see that the warp induces only small deviations about 
dlnPo/dr. In the weakly nonlinear simulation, these deviations are greater and are asym- 
metric about the midplane. This asymmetry is a result of the asymmetric geometry of 
nonlinear warping. In the = slice shown in Figure [2l the "top" of the disk has an 
outward radial pressure gradient (and therefore an inward pressure force), whereas the "bot- 
tom" of the disk has an outward pressure force. Near the disk surface in the warped region, 
the pressure gradients are always larger on the top. These relations are, of course, reversed 
at <p = 7T. As expected, the gradient is weak everywhere in the linear case. On the other 
hand, the gradient in parts of the nonlinear simulation is sharp enough that its scale length 
is ~ 1 in our distance units, ~ r/9 ~ if/2 at the radius shown. Thus, the radial pressure 
gradient is large enough to drive transonic radial motions when ip > 1. 

That these pressure gradients are effective in driving transonic motions is demonstrated 
by the data shown in Figure HI This figure portrays the situation at the same locations 
shown in Figure [3] half an orbit after the beginning of the simulation. In this time, fluid 
elements have traversed a range of 7r in azimuthal angle, sampling the full range of variation 
of the radial pressure gradient shown in Figure [3J In the linear case, there are subsonic fluid 
motions symmetric about zero. In the weakly nonlinear simulation, we see that throughout 
a significant fraction of the disk volume (and a slightly smaller fraction of its mass), the 
initial radial pressure gradients have induced near-sonic radial motions consistent with the 
orientation of the initial pressure gradient. The magnitude of inward and outward radial 
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Fig. 3. — Color contours (see color bar) of dlnP'/dr (Equation [8]) at r = 9 at t = in 
HW-0.5VW (Top) and HW-1.5W (Bottom). No contours are shown where the density falls 
below 0.5% x the maximum. Coordinates are grid coordinates. 
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motion differs due to the asymmetry of the initial radial pressure gradient induced by the 
warp. 

It is these rapid radial motions that mix fluid elements with differing initial orientation 
of angular momentum, ultimately causing the entire disk to align in a single plane. The 
progress of this mixing can be seen in the three panels of Figure [5j There is always a vertical 
gradient in the orientation of the angular momentum due simply to geometry: even when 
all the orbital velocities are parallel, the radius vector from the origin swings with altitude 
across the disk. The magnitude of this purely geometric effect can be seen by the range of 
orientation vertically through the disk at a fixed azimuthal angle 0. However, as the disk 
evolves, the mean orientation at a given <p changes and the range of orientations increases. 
Because there are no external torques on this system, the only way the fluid at a fixed 
location can change the orientation of its angular momentum is by transport of material 
with different angular momentum to that location. This evolution is therefore the signature 
of angular momentum mixing. Moreover, because the range of orientations seen at 1 orbit 
includes directions not present at all at that radius in the initial state, there must be radial 
fluid motions to convey the matter with the new orientation. 

As we have shown, radial pressure gradients associated with disk warping drive transonic 
radial motions; when these streams intersect, their angular momentum is mixed and the warp 
is dissipated. In the traditional approach to this topic, the speed of radial motions is limited 
primarily by a phenomenological isotropic velocity parameterized by a; here, the ultimate 
speed of radial motions is limited by a combination of achieving the maximum speed possible 
for acceleration by a thermal pressure gradient (i.e., Mach number order unity) and the bulk 
viscosity that appears only in regions of strong compression (i.e., shocks). Figure [6] displays 
the rate at which these shocks create entropy. Because strongly nonlinear behavior persists 
for only a small number of orbits in the bulk of the disk, essentially all the entropy production 
is accomplished in the first two orbits. Contrasting different simulations, we find that the 
absolute amount of entropy production rises with both increasing (ip) p and radial extent of 
the warp. 

3.2. Rate of Relaxation Toward the Mean Plane 

The end-result of the relaxation process must be alignment in a single plane because, 
given the constraint of angular momentum conservation and the lack of any external torque, 
that is the only possible long-term equilibrium. Although it is true that in the limit of very 
small amplitude, linear bending waves can propagate indefinitely if there is no dissipative 
process, nonlinear waves — like the ones studied here — create pressure gradients that cause 
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Fig. 4. — Color contours (see color bar) of v r /c s at r = 9 at t = 0.5 in HW-0.5VW (Top) 
and HW-1.5W (Bottom). No contours are shown where the density falls below 0.5%x the 
maximum. Coordinates are grid coordinates. 
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Fig. 5. — Color contours of tan _1 (— L x /L z ) on the shell at r = 14 from simulation HW-1 .5W. 
Three times are shown: the initial condition, 0.5 orbits, and 1 orbit. Only cells where the 
local density is greater than 0.5% of the maximum are illustrated. Coordinates are grid 
coordinates. 
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Fig. 6. — Evolution of volume-integrated entropy for all simulations scaled to its initial value. 
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momentum exchange and, in effect, mix angular momentum. The tilt of the resulting orbital 
plane relative to the equatorial plane must be 6 T {r) = 9 M , where 

Here the local angular momentum is L and the components used in equation [9] are in terms 
of the grid coordinates. To quantify the general level of deviation from this expected equi- 
librium, we define two measures. One is 59 T (r,t), 

5e T (r,t) = e T (r,t)-e M , (10) 

the local angular deviation from the equilibrium state. The other is a mass-weighted spatial 
average of i/j, the normalized radial derivative of the orientation angle 

= l^Mii. 

The time evolution for the latter quantity is displayed in Figured Several clear patterns 
can be seen in this figure. The first is that all the runs exhibit an oscillation in (ip) p . In all 
those runs whose initial state was nonlinear, the amplitude of this oscillation decays rapidly 
at first, but once (4>) p drops below ~ 1, it varies much more slowly for the remainder of 
the simulation; the initially linear simulation oscillates with nearly constant amplitude from 
the start. In other words, $ > 1 is not only an indicator that transonic radial motions will 
be created, it is also a semi- quantitative predictor of the rate at which the warp amplitude 
relaxes: this rate becomes much slower when ijj falls below 1. 

The oscillation is a bending wave induced by the initial warp; in fact, one way of 
describing our results is that an initially nonlinear bending wave rapidly decays to linear 
amplitude. The oppositely directed radial motions created by the warp on the top and 
bottom edges of the disk set up the circulatory flow that characterizes these waves. The 
mode induced is the fundamental (i.e., the radial extent of the disk is half a wavelength) 
because the disk inclination relative to its mean value varies monotonically from one radial 
extreme of the disk to the other. 

It is instructive to model the time- dependent behavior of this oscillation in terms of a 
function f(t) of the form 

fit) = A cos(wt) exp(-st). (12) 

The fitting itself can be done in a variety of ways, however we find that the fit parameters 
(Aq, u, s) depend only weakly on the choice of method. Our method determines the frequency 
of the oscillation uj using the second zero of (if)) p ] we find the parameter s from the magnitude 
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of attenuation at the first minimum of ($) p \ Aq is simply determined by the initial value of 
("0)p. The fit parameters are shown in Tabled with u being replaced by the corresponding 
period, P = 2tt/u, and s in units of inverse orbital periods. 

The most prominent result of this fitting exercise is that over the entire parameter space 
we have studied, relaxation in the nonlinear phase is quite rapid, requiring only a few orbits 
at the mid-point of the disk (the characteristic decay times are ~ 2-5 orbits, rather 
longer than the duration of the relaxation, because the initial values of disk-integrated {i/j) p 
are relatively small, ~l-2, and the rapid decay persists only until (ip) p < 1). The rate 
of relaxation does, however, vary from case to case. If the "severity" of the bend is some 
compound of ip and the radial width over which the disk bends so sharply, the nonlinear 
relaxation rate correlates with "severity" in the way one might expect: more dramatic bends 
relax more rapidly. At fixed Lyy, the fractional decay rate s scales with ip roughly oc ip b , 
with b ~ 1—1.5. When L\y changes, a larger span of warping leads to a larger decay rate; for 
example, the decay rate of HW-2 . 5W is almost exactly twice the rate exhibited by HW-2 . 5N, 
which shares its (ip) p but whose warp stretches only half as far. Conversely, HW-5N and 
HW-2.5W have similar relaxation rates even though their values of (ip) p are a factor of 2 
different; they decay similarly because the case with smaller {ip) p has a warp that is twice 
as wide. A similar pairing exists between HW-2.5N and HW-1.5W. These pairings are also 
evident in the rate of entropy production (see Fig. [6]). 

However, comparing the relaxation rates for HW-1.5W and TW reveals something new 
and unexpected. Following the pattern in which the decay rate scales with if), the rate at 
which the amplitude diminishes in TW is considerably faster than in HW-1 . 5W. However, their 
fractional decay rates s are nearly identical. These two simulations began with identical 
bending rates ip, differing only in sound speed. Therefore, in this regime, the fractional 
decay rate s depends on if), not ip or the sound speed. This is surprising because the primary 
driver of radial motions is unbalanced pressure gradients, so one might have thought that 



Run ID 


A 


P 


s 


HW-0 . 5VW 


0.52 


4.60 


0.12 


HW-1.5W 


1.16 


4.47 


0.22 


HW-2 . 5W 


1.97 


4.33 


0.39 


HW-2 . 5N 


1.07 


4.47 


0.20 


HW-5N 


2.21 


4.33 


0.54 


TW 


3.25 


8.20 


0.21 



Table 2: Quantifying nonlinear warp relaxation. 
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the speed at which fluid elements could mix their angular momenta would always scale with 
the sound speed, yet here the mixing rate appears to have no dependence on the sound speed 
at all. 

That this is possible may have to do with the specifics of our simulations. Typical 
Mach numbers in HW-1.5W during the nonlinear relaxation period are only ~ 0.5-1; at the 
equivalent time in TW, they are ~ 1-3. Both are order unity, but the factor of ~ 3 between 
the characteristic Mach numbers compensates for the factor ~ 3 in the opposite direction 
between their sound speeds. Perhaps if the sound speed were reduced another factor of 3 
(and therefore would require rather greater grid-cell resolution than we can reach at the 
moment), the Mach numbers achieved during the nonlinear relaxation period would not 
be terribly much greater than in TW because it is very difficult for pressure forces alone to 
accelerate motions to more highly supersonic speeds. If so, when ip > 5-10, the mixing speeds 
may reach a limiting value of several times the sound speed. At this point, all we can say is 
that within the range of parameters spanned by our simulations, s ~ O.U(ip/Q.6) b (L w /r)tt 
with b ~ 1-1.5, with no dependence on c s . 

Internal details of the relaxation and bending wave phases can be seen in Figure [HI which 
displays 59r(r, t) and ip(r, t) from simulation HW-1 . 5W. The quantity 59t shows the shape of 
the bending wave more clearly; ip quantifies its degree of nonlinearity. In its initial state, the 
disk bends down at small radii and up at large, reflecting its transition in orientation angle, 
as seen best in S9t- The local normalized bending rate ip, on the other hand, is initially zero 
at both large and small radii and comparatively large in the middle. Because this simulation 
was only weakly nonlinear even in its initial state, ip in most of the disk quickly becomes 
less than unity. As the disk relaxes from its nonlinear warp, it falls into a bending wave 
pattern corresponding to the normal mode whose wavelength is twice the radial extent of 
the disk because this is the mode most closely related to the initial state (once again, seen 
best in 89j). Radially-dependent Fourier analysis demonstrates that these oscillations have 
the same frequency at all radii, confirming that these oscillations constitute a normal mode 
of the disk. Even during the later, predominantly linear phase, both S9t and ip remain 
large at large radii because the disk's inertia decreases sharply toward the edges; regions of 
low inertia have high amplitude because the wave action is conserved along the direction 
of propagation. Wave action conservation alone would lead to nonlinear amplitude at both 
edges, but the nonlinear damping we have already discussed acts on the orbital timescale, 
and is therefore much more rapid at r = 2 than at r = 15. As a result, the wave amplitude 
is relatively small at the inner edge even while it remains nonlinear at the outer edge. 
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Fig. 8. — Spacetime diagram of S9t(t, t) (Top, measured in degrees) and ip(r,t) (Bottom) 
from simulation HW-1 . 5W. 
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3.3. The Angular Momentum Flux 

In order for the disk to arrive at a flat configuration, it must mix angular momentum 
from the different portions of the disk so that all regions have angular momentum with the 
same direction. To describe most clearly how this occurs, it is convenient to define a pair of 
coordinate systems, one spherical and the other Cartesian, both oriented according to the 
mean angular momentum (L) of the disk. To distinguish this Cartesian system from the one 
associated with the grid coordinates, we designate the axes of this one by upper-case letters. 
We will designate the direction of (L) as the polar axis of the spherical system and the Z-axis 
of the Cartesian system. To define the zero-point of azimuthal angle for the spherical system 
and the X-axis of the Cartesian system, we choose a direction that keeps the tilt in the X-Z 
plane. Thus, the dynamics of relaxation can be thought of as the redistribution of Lx of 
differing signs until the local Lx = everywhere. 

Ordinarily, hydrodynamic momentum fluxes are discussed in terms of a Reynolds stress 
tensor R = pv Cg) v . When described in coordinate language, the two velocity vectors are 
nearly always projected onto the same basis. Here, however, we are most concerned with 
the flux of angular, rather than linear, momentum, and the radial motion of the Cartesian 
component Lx- We therefore define the angular momentum flux tensor 

S = J dApv®L, (13) 

— * 

describing v in terms of spherical coordinates, but L in terms of the Cartesian coordinates 
just defined. The integral is taken over a spherical surface at radius r. In this language, 
the element of greatest interest in this tensor is S r x- A natural unit for this tensor is given 
by S (r) = J dArP. We can study the dependence of S r x/S on parameters by looking at 
variations relative to a fiducial case, which we choose to be HW-2 . 5W. 

Figure [9] shows how S r x / So evolves as a function of r and t during the principal relax- 
ation stage in three of our simulations that all have the same warp width: TW, HW-2 . 5W, and 
HW-1.5W. As this figure illustrates, S r x/So is > at all radii for virtually all of the relax- 
ation process, as positive Lx is taken outward from small radii and negative Lx is taken 
inward from large radii. In magnitude, it reaches peak values of ~ 0.5-1, depending on the 
parameters of the simulation. However, these peak values are generally found at the outer 
edge of the disk, not the center of the warp. Near r = r c , the peak value of S r x/S Q ranges 
from ~ 0.15-0.65. It is also clear from this figure that the relaxation takes only a brief time, 
never more than 1-3 orbits at the central radius of the warp. 

The data of Figure |9] can also be used to uncover how the peak magnitude of the stress 
at the center of the warp depends on parameters. Comparing HW-1.5W and HW-2.5W shows 
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that, in rough terms, S r x/So is proportional to the initial warp strength ip. However, the 
fact that the normalized peak stress at r c increases by a factor ~ 4 from HW-1.5W to TW, 
which share the same warp geometry factor ip = d9/d\nr, shows that the peak magnitude 
of the normalized stress is even more strongly dependent upon the sound speed, increasing 
at a rate between linear and quadratic as the sound speed decreases. Put another way, the 
peak unnormalized stress decreases, but only weakly, with decreasing sound speed. That this 
should be so is another reflection of the fact already pointed out in the previous subsection, 
that the Mach number of the radial motions can increase sharply when the sound speed 
decreases. 

We have previously commented on the fact that the exponential decay rate of the warp 
depends on the warp half-width L-\y. This dependence is illustrated in the magnitude of the 
unaligned angular momentum flux. In Figure [TD1 we contrast two simulations with differing 
warp widths, but identical sound speeds and warping rates d6/dlnr, HW-2.5W and HW-2.5N. 
The peak magnitude of S r x/So at r c increases by a factor ~ 2 as the width of the warp 
region doubles, a roughly linear scaling. In other words, the rate at which unaligned angular 
momentum is transported through the disk depends not only on the local warping rate, but 
also on the total extent of the warped region. 

The flux of unaligned angular momentum varies significantly with height within the disk. 
As already shown in Figure HI the radial speeds are greatest on the top and bottom surfaces 
of the disk. Although the density is lower there than in the midplane, Figure [TT1 demonstrates 
that the speed of the radial flows is great enough near the surface to compensate for the 
lower density, making the flux of unaligned angular momentum also greatest near the disk 
surface. In fact, very little unaligned angular momentum moves through the central scale 
height. The direction of the flows is opposite on top and bottom, but the sign of L x is also 
opposite, making S r x have the same sign everywhere. All these patterns are reproduced in 
the other simulations as well. 

It is also worthwhile to examine more closely how the unaligned angular momentum flux 
varies with time during the relaxation. That there should be a non-trivial time-dependence 
is suggested by the fact that the radial pressure gradients are determined by the instan- 
taneous warp rate, but pressure gradients determine instantaneous fluid accelerations, not 
velocities. Reynolds stress, which depends on the velocity, might therefore be expected to 
lag the pressure gradient by some characteristic time. This is indeed the case, as is shown in 
Figure [121 In this figure we plot the time-dependence of both the shell-averaged unaligned 
angular momentum flux at the center of the warp and the ip of the warp shell-averaged at 
the same radius. Time is in orbits at that radius. During the nonlinear relaxation phase, the 
unaligned angular momentum flux S r x noticeably lags the warp rate, while during the later 
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Fig. 9. — The evolution in r and t of the normalized angular momentum flux S r x/So. Three 
simulations with the same Lyy are shown: TW (top panel), HW-2.5W (middle panel) and 
HW-1 . 5W (bottom panel). 
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Fig. 10. — The evolution in r and t of the normalized angular momentum flux S r x/So in 
units of the local pressure and radius. Two simulations sharing the same c s and d9/dha.r, 
but differing in Ljy, are shown: HW-2.5W (top panel) and HW-2.5N (bottom panel). 
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linear bending wave phase, the lag grows still longer. For fixed sound speed, the magnitude 
of this lag depends primarily on the width of the warped region, not the warp rate. For 
example, as shown in Figure IT2"} the lags in the two runs with Lyy = 2 (HW-5N and HW-2 . 5N) 
are almost identical, both ~ 0.6 orbits, while the lags in the two runs with L\y = 4 and the 
standard sound speed (HW-2.5W and HW-1.5W) are likewise almost identical at ~ 0.8 orbits. 
On the other hand, TW, which also has L\y = 4, but a sound speed smaller by a factor of 2.5, 
shows a lag ~ 1.4. Thus, the characteristic acceleration coherence time appears to be very 
roughly given by ~ l.§{L w /c s ) 1 / 2 ~ A[L W / (cstt)} 1 ' 2 . Moreover, this same lag between warp 
rate and angular momentum flux also explains the overshoot that creates a persistent linear 
bending wave after the nonlinear warp rate has been eliminated. 

In those cases in which the time- dependence of S r x approximately follows the time- 
depence of ip at a well-defined lag r, it is of further interest to ask how S r x(r,t)/ S depends 
on ip(r,t — t). We have done so, but find it is rather complex. At any given radius, there 
is generally a clear relation between the two in which S r x/So increases with increasing ijj, 
but it is typically strongly nonlinear. In addition, the character of that nonlinearity varies 
substantially with radius within a single simulation, as well as from case to case. Thus, 
these data do not reveal any consistent functional relation between S r x and if), only that in 
a rough qualitative sense one rises as the other does. 



3.4. Comparison with Analytic Formalisms 



As we have already remarked, the parameters of our study do not match well with the 
regimes in which any previous analytic work is valid. Unfortunately, these previous analytic 
efforts are the only ones with which to compare our work, so they are the only options for 
any attempt to provide some context. In this section we attempt to situate our work relative 
to that context. 

In past treatments of warp dynamics, the radial mi xing of angular momentum by warps 
has most often been modeled by analogy with di ffusion (Pringlel 19921; Ogilvie 1999 ). even in 
cases in which a is not large compared to H/r ( Lodato fc Price! 201Clh . Two reasons make 
it worthwhile to contrast our results with those of the diffusion picture, even though our 
simulations are formally in what is us ually called th e "bending wave" regime. One is that 
the weakly nonlinear diffusion theory of lOgilvid ( 119991 ) was applied to this case, but was found 
to fail for the specific parameters we study here; comparison with that theory's assumptions 
may therefore be instructive. The second is that the physics underlying the diffusion model 
is in fact not so different from what we have studied. In both pictures, warps create radial 
pressure gradients that drive radial flows, and these radial flows mix angular momentum 
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radially. Where they differ is primarily the mechanism that limits the velocities of these 
flows and secondarily the nature of the microscopic diffusion that completes the angular 
momentum mixing. In the diffusion model, a phenomenological shear viscosity restricts the 
magnitude of radial motions because their speeds vary with height from the disk plane; in 
our study, the limit is imposed by a combination of the finite extent of the radial pressure 
gradient and shocks. In the diffusion model, that same phenomenological shear viscosity 
directly mixes momentum; in our case, mixing is achieved by a combination of the localized 
artificial bulk viscosity mediating shocks and numerical diffusion when velocity gradients are 
steep even at the grid scale. 

As is implied by results we have already discussed, what we observe resembles diffusion 
in certain ways, but not in all. Like diffusion, at a single radius at a given time there are 
fluid elements flowing both in and out, with any net fluxes attributable to the remainder 
after these flows nearly cancel. Similarly like diffusion, the end result is to mix conserved 
quantities, in this case angular momentum. However, there are also respects in which this 
process differs from diffusion. For example, there is comparatively large-scale organization 
of the regions flowing inward and outward, and the net displacements of fluid elements are 
not particularly small with respect to the gradient scale. Nonetheless, for the purposes of 
mathematical modeling, all that really counts is whether the flux of angular momentum 
can be described as a diffusion coefficient times the gradient in angular momentum. More 
precisely, accounting for the net inward motio n of material due solely to the diffusion of the 



unaligned angular momentum (iPringld Il992l ). a diffusion model predicts that the rate at 



which unaligned angular momentum passes through the shell at radius r is 

„ 2o / . sine de \ de 
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where v 2 is the diffusion parameter. With our simulation data, we can test whether the flux 
is related to the bending rate in this manner. In particular, the key property that defines a 
diffusion-like model is whether the flux at a given time and location can be described as the 
product of a diffusion coefficient (possibly dependent upon the warp) and the warp rate, both 
evaluated at the same time and location as the flux in question. For our simulations, the 
maximum departure angle from the mean angular momentum orientation is ~ 0.3 radians, 
so the second term inside the parentheses in equation [TH is small. 

Viewing the results reported in the previous subsection from the point of view of this 
question, several additional contrasts with the diffusion picture appear. One is that the 
magnitude of the unaligned angular momentum flux S r x responds not just to the local 
orientation gradient d6/d\nr, but also to the width Lw over which the disk is warped. 
Another is that there appears to be a significant delay between the time-dependence of 
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S r x and the time-dependence of d9/dhxr, with this delay scaling in rough proportion to 
[Lw I \c s Vl)] 1 / 2 . For the parameter range explored here, the magnitude of this delay is ~ 
1 orbit, comparable to th e nonlin e ar wa rp relaxation time. By contrast, in the asymptotic 



expansion at the heart of lOgilvid ( 119991 ). it is assumed that the shape of the disk changes 
only on timescales ~ (H/r)~ 2 x longer than the orbital time. Still another is that although 
S r x consistently increases with dO/dlnr when the warp rate is evaluated at an appropriately- 
chosen earlier time, the data do not suggest any consistent function al relation betw een the 



two. It is conceivable that the inability of the nonlinear model of lOgilvid (119991 ) to find 
a solution for «2 in cases whose parameters are similar to ours (zero anomalous viscosity, 
epicyclic frequency smaller than the orbital frequency) is related to these disparities. 

On the other hand, as we have already remarked, the elimination of any anomalous 
viscosity from our model may be taken as an indication that the nonlinear warp problem 
we have treated resembles linear bending waves more than linear warp diffusion. There 
are both significant para llels and significant departures between linear bending wave theory 



( jLubow fc Ogilvid 120001 ) and what we observe here. That the peak stress increases approx- 
imately in proportion to the bending rate would be natural for a linear theory. Similarly, 
a delay between the bending rate and the torques it induces lies at the heart of bending 
wave theory. However, for linear bending waves, the delay timescale defines the period of 
the oscillations, not the time required for the initiation of damping. In fact, when a = as 
in our simulations, linear bending wave theory predicts that there is no damping at all. By 
contrast, we find that when the bend is nonlinear, it is damped by exactly those Reynolds 
stresses responsible for wave propagation when the bend is weak. Further contrasts can be 
seen in two other facts: that the magnitude of these Reynolds stresses increases in rough 
proportion to L\y, effectively the wavelength of the bending wave; and that the normalized 
Reynolds stress increases rapidly with decreasing sound speed. Such nonlinear couplings are, 
of course, entirely absent in linear theory. 

The concentration of S r x near the disk surfaces also underlines the importance of treat- 
ing properly stresses acting on the vertical shear flow. Here we have taken an extreme posi- 
tion, assuming zero anomalous viscosity. Most treatments assume that there is an isotropic 
anomalous viscosity, one that responds to r — 9 shear as strongly as to r — shear. Real 
internal stresses due to MHD turbulence might act in a different way, both in terms of their 
relationship to shear and in terms of their directional dependence. 
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4. Consequences and Conclusions 

We have studied a series of simulations to assess the evolution of strong warps in inviscid 
disks. We find that the evolution of these disks is governed by the parameter tp, the radial 
extent of the warp L w , and the sound speed c s . A finite thickness flat disk with finite 
radial extent can achieve hydrostatic equilibrium by pressure gradients with both radial and 
vertical components: the radial components are balanced by (usually small) departures from 
Keplerian rotation, while the vertical components are balanced by the vertical component of 
gravity. However, a warped disk cannot be in such an equilibrium; it must contain unbalanced 
radial pressure gradients. As a result, purely acoustic effects stemming from these pressure 
gradients force radial motions. When > 1, these motions are transonic, mix angular 
momentum radially, and rapidly relax the warp. 

In our simulations, conducted with neither any explicit viscosity nor the sort of MHD 
stresses that convey angular momentum in a conventional flat disk, the speeds of the radial 
motions depend on the overall strength of the warp, a combination of (ip) p and the radial 
extent over which the bending takes place. In the parameter regime we studied, the non- 
linear exponential damping rate s ~ 0.14(ip /O.Q) b (Lw /r)Q with b ~ 1-1.5, and is almost 
independent of c s . When we changed the sound speed at fixed ip, the Mach number of the 
radial motions changed in the opposite direction, leaving the radial speed largely unaltered. 
It is possible that the Mach number saturates for still larger values of ■?/>, so that s declines 
with decreasing c s , but in the parameter range covered by our simulations no such trend 
was apparent. Even though the fluid motions are only transonic, in disks like the ones we 
have studied, where H/r is not <C 1, it is possible for streams to cross the entire span of the 
warp in roughly an orbital period because (r/c s )(Q/2ir) ~ (l/2ir)(H/r). As a result, in our 
simulations relaxation from a state of strongly nonlinear warping to linear behavior required 
only a few mid-warp orbits. 

Relaxation of warps in a purely hydrodynamic disk requires fluid motions to mix angular 
momentum from regions initially having different orientations. The flux of unaligned mo- 
mentum in the radial direction, a quantity we dub S r x, therefore becomes the key quantity 
governing the relaxation. Scaling this quantity in units of the product of the local pressure 
and radius, we find that S r x increases with warp rate, but also with warp width Lw] in other 
words, it is controlled in part by global conditions, not local. Moreover, as might be expected 
when the fluid motions carrying angular momentum arise from the acceleration in pressure 
gradients directly associated with the warp, S r x varies in a way that is related to the warp 
rate, but delayed by a lag that scales oc (Lw / c s ) l l 2 . For the parameters of our simulations, 
these lags were comparable to the decay time; in cooler, thinner disks, this scaling suggests 
that they become even longer. 
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These results suggest that decay of nonlinear warps conforms to neither a conventional 
diffusion nor a bending wave formalism. The process depends on a net divergence in un- 
aligned angular momentum flux, as in a diffusion picture, but this flux is determined both 
by the history of the warp (because of the significant delay between changes in the warp and 
the creation of angular momentum fluxes) and global disk properties. On the other hand, 
bending wave dynamics do not entirely suffice because the same Reynolds stresses that lead 
to bending wave propagation in the linear regime create rapid warp relaxation, i.e., wave 
damping, when they become nonlinear. Moreover, no anomalous viscosity is required to 
control the speed of the radial motions; purely hydrodynamic effects, such as the inability 
of thermal pressure gradients to accelerate fluid to more than a few times the sound speed, 
and the creation of weak shocks, suffice. 

The simulations presented here are by no means the first computational study performed 
of warped disks. There has been previous work ( jNelson fc Papaloizoulll999l ; lLodato fc Pringle 
20071 ; lLodato fc Pricell2010l ) using SPH simulations to study this phenomenon. However, our 
simulations are both the first to use a global grid- based treatment and the firs t to focus on 
hydrodynamics without any anomalous viscosity. iNelson fc Papaloizoul (119991 ) studied the 
propagation of bending waves excited by pulses of varying amplitudes placed in the mid- 
dle of a disk whose imposed anomalous viscosity was small compared to the disk thickness 
(a < H/r). Like them, we found strong damping of nonlinear amplitude disturbances by 
radial motions, but they did not study how the relaxation rate o r its associated angular 
momentum fluxes dep ended on parameters. More recent SPH work (lLodato &: Pringldl2007l ; 
Lodato fc Pried 120101 ) has focused on the "diffusive" regime, in which a > H/r; even the 
simulations they labeled as "inviscid" assumed an«~ 2H/r. They are therefore not directly 
comparable to our work. In addition, like the earlier simulational efforts, they did relatively 
little in the way of parameter exploration or investigation of the detailed time-dependence 
of warp relaxation, so we cannot compare our results in those regards to theirs. Nonethe- 
less, it is interesting to note that lLodato fc Pried (120101 ) speculated that when a is small 
and d6/d In r comparatively large (i.e., in our circumstances), a diffusion model might be 
inappropriate. 

The absence of ordinary internal stresses (i.e., those produced by MHD turbulence or 
often modeled by an anomalous viscosity) might be considered a limitation of this work. 
We chose to eliminate them in order to clarify which processes account for which effects. 
In future work on this subject we will include explicit MHD stresses in order to see what 
new effects they introduce. In particular, it will be especially important to see the degree to 
which MHD turbulence restricts the growth of the transonic radial motions crucial to warp 
relaxation, and what relation that degree of restriction has to the stresses responsible for 
angular momentum transport in flat disks. As we have shown, the concentration of unaligned 
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angular momentum flux to the surfaces of the disk is so strong that careful treatment of any 
stresses induced by this shear could be quite important. It is possible, however, that the 
relevant MHD stress (the r — 9 component) will remain a small effect. In most simu l ations 
of MHD turbulence stirred by the magnetorotational instability from IStone et al.l (119961 ) 
onward, Bg is the smallest of the three field components, only ~ O.IB^. Consequently, even 
at distances one or two scale heights off the midplane, where the shear is strongest, the MHD 
stress restraining the radial motions may be weak. 
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A. Numerical Tests 



The rotational invariance of the Euler equations implies that the simulation results 
should be invariant under a rigid rotation. Although this is true physically, this invariance 
will not, in general, be respected numerically. In special cases, when the velocity is parallel to 
a symmetry direction of the flow and this direction is along a coordinate axis, fluid transport 
creates very little numerical dissipation. A statistically axisymmetric orbiting disk treated 
in either cylindrical polar or spherical coordinates whose equatorial plane coincides with the 
disk plane is an example of this favorable special case. When that same disk is inclined 
to the equatorial plane, however, the numerical dissipation can be considerably enhanced. 
Unfortunately, when a disk's inclination changes from place to place, it is impossible to avoid 
having the dominant velocity oblique to the coordinate system somewhere, yet understanding 
the physics associated with a warped disk relies on being able to separate physical effects 
from numerical ones. We therefore took special pains to choose a combination of coordinate 
system, inclination, and grid resolution that would keep numerical dissipation small enough 
to be negligible. 

Our procedure was to consider a series of disks with a constant, but non-zero, inclination 
angle in both cylindrical and spherical coordinates at both low and moderate resolution, 
roughly 4 ZPH and 8 ZPH respectively (ZPH is zones per scale height). Both grids employed 
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isotropic cells. These flat, but tilted, disks had the same radial extent and thickness as 
the production run initial conditions. Each simulation was run for 10 orbits at r = r c . To 
quantify numerical losses, we computed two measures: the fractional change in the magnitude 
of the fluid's total angular momentum A|L|/|L| and the fractional change in the density- 
weighted tilt 9t- The latter quantity is defined as 

6 T (t) = tan" 1 ( ~ < L *> \ . ( A1 ) 



< L z > 

Here the components of the angular momentum L are taken in the grid coordinate system. 
The results from this study are shown in Table [3j Note that the net change in the magnitude 
of the angular momentum is always negative, and that any change in mean inclination is 
always in the sense of bringing the disk closer to the coordinate system's equatorial plane. 

Comparing cylindrical and spherical tests at the same 9t, we find that spherical coor- 
dinates are clearly better according to both measures when 6t < 45°. At 6? = 45°, the 
verdict is mixed: spherical coordinates are better by the AL/L criterion, cylindrical by the 
A6t/0t criterion. At larger angles, both coordinate systems do poorly, and to a similar 
degree. The results for low inclination angles, however, give us confidence that the study 
of warped disks in which the inclination angle doesn't exceed the threshold of 45° can be 
conducted with confidence, provided the resolution is at least 8 ZPH. For our production 
runs, we used exclusively spherical coordinates due to their superior performance for the 
inclination angles of interest. In particular we note that our warped simulations are most 
closely comparable to the oblique disk simulation S30M, for which the fractional losses in 
both the angular momentum and inclination angle are less than or approximately 1%. 

There exists much work in the literature discussing criteria for the proper resolution 
of an accretion disk simulation. Although satisfying these constraints will be particularly 
important when internal stresses studied because they are caused by MHD turbulence, for 
pure hydrodynamics the resolution constraints are much weaker. Nonetheless, we tested 
the resolution-dependence of our results. We conducted two simulations identical save for 
their resolution, which corresponded to roughly eight and sixteen zones per vertical scale 
height. Snapshots of the evolution with time of their inclination profiles, 9t(t, t) are shown 
in Figure [131 where 

6 T (r lt ) = tan' 1 f ~ >6 A , (A2) 



and the notation X(r) =< X > d ^ is used to denote a shell average of a quantity X. Overall, 
we find excellent agreement, and in particular find almost indistinguishable behavior within 
the body of the disk (2 < r < 15). The comparison shown here is limited to 2.5 orbits at 
r = r c , but, as we will show, it is during this period that the most violent realignment of 
this disk occurs. This is therefore the most important time to examine. 
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Phi/x 

Fig. 11. — The unaligned angular momentum flux S r x in HW-5N at the middle of the warp 
(r = 9) in the middle of the nonlinear relaxation phase (t = 0.5 orbits) as a function of 
azimuthal angle and polar angle 9 from the axis. 



Run ID 


-AL/L 


-A6j< 1 6x 


C15L 


1.45% 


5.46% 


C15M 


0.63% 


1.83% 


C30L 


3.24% 


6.76% 


C45L 


11.08% 


7.44% 


C45M 


3.51% 


1.99% 


C60L 


23.63% 


9.20% 


S15L 


0.54% 


2.58% 


S15M 


0.41% 


1.07% 


S30L 


1.49% 


3.20% 


S30M 


0.71% 


1.12% 


S45L 


8.56% 


7.28% 


S45M 


3.29% 


2.46% 


S60L 


24.61% 


12.74% 



Table 3: Results from a series of runs to test the numerical importance of disk inclination. 
Simulations are denoted by the geometry used: (C)ylindrical or (S)pherical; the angle (in 
degrees) of the oblique inclination; and whether the resolution was (L)ow or (M)oderate, 
using 4 or 8 zones for each vertical scaleheight. 
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Fig. 12. — Time-dependence of the shell-averaged t[j (solid curve) and AS rX /S a t the mid- 
dle of the warp region in (from left to right and then top to bottom): HW-HW5N, HW-2.5W, 
HW-2.5N, HW-1.5W, and and TW. With this arrangement, the two narrow warp width simula- 
tions are on the left, while the wide warp width simulations are on the right and the bottom. 
The normalized stress is multiplied by 4 so that its time-dependence can easily be compared 
with that of ■0. Time is in units of orbits at r = 9 in the first two cases, at r = 8.3 in the 
third. 



-36 - 




Fig. 13. — Resolution test of a warped disk. Shown are the radial inclination profiles at 
various times for two simulations differing only in their resolution (8 ZPH: dashed, 16 ZPH: 
solid). 
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